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Abstract We consider here the effects of inertia on 
the instability of a flat liquid him under the effects of 
capillary and intermolecular forces (van der Waals inter¬ 
action). Firstly, we perform the linear stability analysis 
within the long wave approximation, which shows that 
the inclusion of inertia does not produce new regions 
of instability other than the one previously known from 
the usual lubrication case. The wavelength. Am, corre¬ 
sponding to he maximum growth, ujm, d.nd the critical 
(marginal) wavelength do not change at all. The most 
affected feature of the instability under an increase of the 
Laplace number is the noticeable decrease of the growth 
rates of the unstable modes. In order to put in evi¬ 
dence the effects of the bidimensional aspects of the flow 
(neglected in the long wave approximation), we also cal¬ 
culate the dispersion relation of the instability from the 
linearized version of the complete Navier-Stokes (N-S) 
equation. Unlike the long wave approximation, the bidi¬ 
mensional model shows that Am can vary significantly 
with inertia when the aspect ratio of the him is not suf¬ 
ficiently small. We also perform numerical simulations 
of the nonlinear N-S equations and analyze to which ex¬ 
tent the linear predictions can be applied depending on 
both the amount of inertia involved and the aspect ratio 
of the him. 


1 Introduction 

The stability of thin films on substrates has been for a long time a basic subject of 
research, not only because of the numerous technological applications, including coat¬ 
ings, adhesives, lubricants, and dielectric layers, but also because of their fundamental 
interest (Craster and Matar, 2009; Eggers, 1997; Oron et ah, 1997). The dewetting of 
thin liquid films is the process of destabilization of such films which leads to the for¬ 
mation of drops. It is generally observed when the supported liquid film is placed on 
a substrate under partial wetting conditions, and subject to destabilizing intermolecu¬ 
lar forces. For a homogeneous isotropic liquid on a uniform solid substrate, two main 
dewetting processes are known: (i) the nucleation of holes at defects or dust particles, 
and (ii) the amplification of perturbations at the free surface (e.g., capillary waves) un¬ 
der the destabilizing effect of long-range intermolecular forces in the so-called spinodal 
dewetting (Thiele, 2003; Thiele et ah, 2001, 1998). Although the distinction between 
these two dewetting processes is well established in the literature, there is still a lot of 
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debate about which of these mechanisms is actually observed in a given experiment. 

In this context, lubrication approximations to the full Navier-Stokes equations have 
shown to be extremely useful for investigating the dynamics and instability of thin liquid 
films on substrates, including the motion and instabilities of their contact lines (Oron 
et ah, 1997). The theoretical treatment of the coating problem is greatly simplified if the 
film is so thin that the lubrication approximation can be employed. When this modeling 
is valid, it is possible to determine the velocity field of the liquid as a function of the 
film thickness, and the problem reduces to the solution of a nonlinear evolution equation 
for the thickness profile of the film. To leading order, at low speeds, the dynamics is 
controlled by a balance among capillarity, viscosity, and intermolecular forces, without 
inertia playing a role. This approach has achieved considerable success in dealing with 
the solution of this class of problems (Colinet et ah, 2007). 

However, in some applications such as the dewetting of nano-scale thin metallic films 
on hydrophobic substrates, the effects related to inertia and the shortcomings of the 
lubrication approximation assumptions (requiring small slopes and consequently small 
contact angles) appear to have a crucial influence on the dynamics and morphology of the 
film (Gonzalez et ah, 2013). Experimental studies of unstable thin films coating solids 
have shown significant differences in the patterns that develop when fluid instabilities 
lead to the formation of growing dry regions on the solid. The effects of inertia on 
the instability have been studied previously in other problems, for example for a film 
flowing down an incline (Lopez et ah, 1997), the breakup of a liquid filament sitting 
on a substrate (Ubal et ah, 2014), and several other configurations (Hocking and Davis, 
2002). However, these problems do not include explicitly the effects of the intermolecular 
interaction between the molecules of the liquid and those of the solid. Here, we consider 
it by using integrated Lennard-Jones forces, which lead to the disjoining pressure that 
entails the power dependence on the fluid thickness (Kargupta et ah, 2004). In the 
present context, the occurrence and nature of both inertia and bidimensional effects in 
the liquid film on the solid substrate is of interest, not only for fundamental research, 
but also for technological applications. 

The solutions of problems under the lubrication approximation is usually limited to 
speeds low enough to give small capillary and Reynolds numbers. The extension of 
the theory to higher speeds introduces inertia into the problem, and, even in the case 
of thin films, the analysis may become much more difficult. The great simplification 
previously found by the application of the lubrication theory no longer exists; instead, 
the system is governed by the coupling of a nonlinear partial differential equation for the 
velocity field, and an evolution equation for the thickness profile. It is possible, however, 
to find a class of problems in which inertial effects can be assessed within the long 
wave framework. In this work we are concerned with the instability of a flat liquid film 
extended over a solid plane, and subject to intermolecular forces between the liquid and 
the solid substrate. Then, the film evolution is studied by considering viscous, surface 
tension, and intermolecular forces, with special emphasis on the effect of inertia in the 
development of the instability. 
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2 Intermolecular forces in the hydrodynamic description 

We consider a thin liquid film of thickness /iq, which spans infinitely in the x-direction 
(the system is invariant in the y-direction, i.e. plane flow conditions prevail), and rests on 
a solid plane at z = 0 (see Fig. la). Here, we will consider the instability of this initially 
flat film under the action of surface tension and intermolecular forces, acting both at the 
free surface of the film with instantaneous thickness h{x,t). Thus, the hydrodynamic 
evolution is governed by the Navier-Stokes equation and the incompressibility condition. 


p {dt'v + V • Vv) = —Vp + pAv, V • V = 0 


( 1 ) 


where p is the liquid density, p, its viscosity, p the pressure and v = (rt, w) the velocity 
field. At the substrate (z = 0), we apply the no-slip and non-penetration conditions. At 
the free surface, z = h{x,t), we have the usual kinematic condition and normal stress 
equilibrium given by 

p = -U-^C ( 2 ) 

where 7 is the surface tension, C the curvature of the surface, and 


U{h) 


Kf{h) = K 




( 3 ) 


is the disjoining pressure. Here, k is a constant with units of pressure (related to the 
Hamaker constant of the system), the exponents satisfy n > m > 0, and h* is the equilib¬ 
rium thickness (of the order of some nanometers). This surface force is a consequence of 
the interaction among the molecules in the three phases present in the problem, namely 
the liquid of the film, the solid substrate and the surrounding gas. Note that at equi¬ 
librium, i.e. when h = ho = const., the film has a uniform pressure po = ~n(^o) > Oj 
since n(/io) < 0 for /iq > /i*. 


3 Long wave approximation 


In this approximation it is assumed that the film thickness, /iq, is much smaller than 
the characteristic horizontal length of the problem. Since the film extends to infinity, we 
assume that there exists a typical length associated with the wavelength of the perturba¬ 
tions, namely i. The definition of i will be made more precise below. Subsequently, for 
e = /io/£ <C 1, we can simplify Eq. (1) under the long wave approximation assumptions 
retaining inertial terms in the form 


dp d‘^u 
dx ^ dz"^ 
dp 
dz 

du dw 
dx dz 


( du 


du 


du 


= P\ nr 


\dt 


dx 


dz 


-tA = 0 


= 0 . 


( 4 ) 

( 5 ) 

( 6 ) 
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Figure 1: (a) Schematic diagram of the problem, (b) Parameter A as given by Eq. (24) 
as a function of the ratio between the equilibrium thickness, /i*, and the film 
thickness, ho, for two pairs of the exponents {n,m). The vertical dotted lines 
correspond to go = 0, i.e. /i* = ho 


The boundary conditions for these equations are as follows. At 2 = 0, we impose no 
penetration and no slip at the substrate. 


re = 0, u = 0. 

At the liquid-gas interface (z = h), we have zero shear stress, 


as well as the kinematic condtion. 


^ = 0 

dz 


dh dh 


and the Laplace relation for the capillary pressure 


d'^h 


pW = - 7 ^ - 


( 7 ) 


( 8 ) 


( 9 ) 


( 10 ) 


From Eq. (5) we see that the pressure, p, is z-independent, and then p is only a function 
of h, p = p{h). Thus, we have that the x-derivative of p in Eq. (4) is given by 


dp d^h f,^,^dh 


( 11 ) 


The continuity equation, Eq. (6), is conveniently satisfied by introdneing the stream 
function '0(x, t) defined by 


u = 


dip 

Ih' 


w = — 


dip 
dx ' 


( 12 ) 
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Therefore, Eqs. (4) and (9) in terms of il) are given by 


^ dz^ dx^ ^ dx~^ ^ \dzdt dz dxdz dx dz'^ ) ’ 


^ d'il;{x,h,t) ^ 
dt dx 

The boundary conditions, given by Eqs. (7) and (8), in terms of ip are: 


V’L=o = 0> 


dip 

dz 


= 0 , 

z=0 


d'^ip 


= 0 . 

z=h 


(13) 

(14) 


(15) 


3.1 Linear stability analysis within long wave approximation 

The equilibrium state is given by /i = /iq, and for small-amplitude perturbations, the 
height and stream function can be written in the form 


h = hQ(l + , V’ = AiPi{z)e‘^^+^'^^, 


(16) 


where A is the small amplitude of the perturbation, and unstable (stable) modes corre¬ 
spond to a; > 0 (w < 0). By replacing Eq. (16) into Eqs. (13) and (14), and retaining 
terms up to order one in e, we have 


(fipi 

dz^ 

uj -|- ikipi{hQ) = 0 , 


/i ^ = ijhok^ — iKhof'{hQ)k + poj 


dipi 

dz 


with the boundary conditions 


II n 

-/-A-o-o. ^ 


= 0 , 


d^ipi 


z=0 


dz'^ 


= 0 . 


2 = 1 


(17) 

(18) 


(19) 


Now, we define the horizontal length scale, i, by choosing Kf'{ho) = 7 /^^, so that it 
turns out 

where 

90 = =-n (I) +”'Cip) ^ (21) 

Since n > m and ho > h^:, we have go > 0. Note that we are here including in I all the 
effects related with the magnitude of the intermolecular forces given by k. In fact, this 
constant is usually related in the literature with the contact angle, 0, which appears at 
the contact regions formed when the film thins up to /i*, and characterizes the partial 
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wetting of the substrate. It is found that the following simple relationship holds (Oron 
et ah, 1997; Schwartz and Eley, 1998; Diez and Kondic, 2007) 


K = 


7(1 — cos 9) 
M/i* 


( 22 ) 


where M = {n — m)/{{n — l)(m — 1)). Thus, the characteristic length scale becomes 


Mhoh^ 


(1 - cos 9) go' 


(23) 


so that this length includes all the parameters determining the problem, except for 7 
and g which yield the time scale (see Eq. (26) below). In Fig. lb, we show how the 
dimensionless combination 


A = Vl 




(24) 


depends on the ratio h^/ho for two fixed values of the exponents pair (n,m). Inter¬ 
estingly, very small values of h* as well as /i* close to ho {m /n) yield very large 

values oi I j ho for given contact angle, 9 < 7 r/ 2 . 

Consequently, a convenient non-dimensional version of the problem for the long wave 
approximation is given by the following scaling 

A-=f. Z=A = K = (25) 


where 

ni 

7 

is the time scale. Under these definitions, Eqs. (17) and (18) become 


(26) 




(27) 


Q + iK^i{l) = 0, 

(28) 

where 


(29) 

with 

La* = La 

(30) 

and 

T - 

,,2 

(31) 


being the Laplace number. The latter dimensionless number considers the effects of all 
the forces playing a role in the flow, namely, inertial (characterized by p), viscous (char¬ 
acterized by p), surface (characterized by 7 ) and intermolecular forces (characterized by 
£). 
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The solution of Eq. (27) has the form 

>1.. = iK{K-^ - i) <lZ + sm(9-,Z)secZ-tan., _ 

which allows one to obtain the dispersion relation from Eq. (28) as 

0 q — tan q 

A2 (772 _ 1) = ^3 • 

In the limit g —?• 0, this expression tends to the purely viscous solution, 

O K\l-K^) 


(32) 


(33) 


(34) 


which is obtained in the inertialess case (Diez and Kondic, 2007) for La = 0. Note 
that the dimensionless critical (marginal) wavenumber is equal to unity for the viscous 
case, i.e. Kc = 1, so that Eq. (34) shows instability for A < 1. This is because the 
choice of the in-plane characteristic length, i, the inverse of the dimensional critical 
wavenumber (Nguyen et ah, 2012). 

By dividing Eq. (33) by and using Eq. (29), we may dehne the parameter r as 


1 


r = 


La* A2(A2 - 1 )’ 

and then, the possible values of q for given K, are given by the roots of 

tan q — q 


r = 


(35) 


(36) 


In what follows, we will consider only real values of K. Thus, the allowed values of r 
are r < rmax = —4:/La* for A < 1, and r > 0 for A > 1 (see Eq. (35) and Eig. 2a). 
In the region A < 1 and r < Tmax there exist two different values of A for a given r, 
and so they share the same growth rate, D. At r = rmax we have A = = l/\/2- 

Instead, in the region A > 1 and r > 0, each mode A has a unique and different r, and 
consequently, D. 

In order to analyse the possible values of D in each region, it is convenient to introduce 
the notation 

(37) 


q = 


so that the complex growth rate is 


La* 


(38) 


Eor Llr = 5R(D) > 0 (< 0) we have unstable (stable) modes, and for Dj = ^(D) 7 ^ 0 we 
have spatially oscillating modes. Therefore, we consider the imaginary and real parts of 
Eq. (36), which read as 


5R(r) = F{qo, If) = [—<1> cos4(/? + sin(2go cos (^) cos 5(/? + sinh(2qo sin(^) sin5(/?]/A, (39) 
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K 


(a) La*r vs. K 



Figure 2: (a) Relationship between La*r and K as given by Eq. (35). The dashed line 
is r = —4:/La*, (b) Possible values of r as a function of go- The blue lines 
correspond to (/? = 0, the red one to (/? = 7 r /2 and the black one to 7 ^ const, 
(see Fig. 3a). 


^(r) = G(g 0 ) sin — sin(2go cos ip) sin bp + sinh(2go sin p) cos bp] /A, (40) 

where 

$ = 2go cos (goe“*‘^) cos (goe*‘^) , (41) 

A = gQ(cos[2go cos (/?] + cosh[2go sin y?]). (42) 

Since r is real, the solntions of Eq. (36) must have ^(r) = 0. Two trivial roots of this 
function are (/J = 0 and \p\ = 7 r/ 2 . However, it is possible to find roots also along a curve 
in the {qo, p) plane given implicitly by the function G{qQ, p) = d (see Fig. 3a). 

For \p\ = 7 r /2 we find unstable real modes with growth rates given by (see Eq. (38)) 

H = ^ > 0, (43) 

LjOl 

where go is now given by the implicit relation 

r{±i go) = F{±i go, ± 7 r/ 2 ) = ( 44 ) 

% 

The function r(±fgo) is plotted in with red lines in Fig. 2b. Since r < 0 for all go, this 
branch corresponds to AT < 1. 

Instead, for (^ = 0, we obtain stable real modes whose growth rates are given by (see 
Eq. (38)) 

VL = VLr = —— < 0, (45) 

La* 

where go is obtained through the implicit relation (see blue lines in Fig. 2b ) 

r(go) = F(go, 0 ) = --- ^V'^° - (^ 6 ) 

Qo 
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% 



1 / La* 


(a) 


(b) 


Figure 3: (a) Curves in the {qo, (p) plane along which 3'(r) = 0. (b) Maximum growth rate 
as a function of La*~^. The horizontal dashed line corresponds to ^vis,max = 
1 / 12 . 


The function r(g'o) is plotted in with blue lines in Fig. 2b. Since r changes sign at 
qq = 7r/2, the upper branch corresponds to AT > 1, while the lower one to K < 1. 
Moreover, these branches are related to monotonically damped modes. 

The implicit relation G{qo,(p) = 0 (plotted as the black curve in Fig 3a) allows to 
obtain r(goe*‘^) as a function of qo (see black curve in Fig. 2b). This branch appears as 
a bifurcation point of the the upper branch ip = 0, with coordinates B = (r;,, qo^j,) = 
(1.1127,0.5367). Since \ip\ 7 ^ 0,±7r/2 we have complex values of the growth rate, Q, 
as determined by Eq. (38). Moreover, since \ip\ < 7r/2, is always negative and it 
corresponds to oscillating stable (damped) modes. Besides, it turns out that r > 0, so 
that these modes belong the region AT > 1. 

As a result, only the branch \ip\ = 7r/2 includes unstable modes, which are in the 
region k < 1 and r < Vmax of Fig. 2a. The mode with maximum (real) growth rate, 
^maxi is given by r = rmax < 0 for a given La* and is located at the intersection with 
the line \ip\ = 7r/2 in Fig. 2b. In fact, for given La* we solve 

tanh qo^max q0,max 4 (AV) 

'dmax 

for qo^rnax and obtain Ll^ax = ~Qo max /result is shown in Fig. 3b, where it is 
observed how Clmax tends to the viscous value, namely Llms,max = 1/12 (see Eq.(34)), as 
La* —>-0. It is also shown that the behavior for large La* corresponds to a decreasing 
growth rate as a power law with exponent close to 0.42. Similar decreasing trends of the 
growth rates due to inertial effects have also been found in other problems (Oron et ah, 
1997; Ubal et al, 2014). 

Eigure 3b also shows that the line r = Vmax is also intersected by the (/? = 0 line. 
Since it corresponds to monotonically damped perturbations in the region 0 < /c < 1, 
this implies that the maximum damping for the stable mode occurs at the same k than 
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the unstable modes in the \ip\ = 7r/2. 

Note that unstable monotonically growing modes are only possible for A: < 1, so that 
neither the critical wavelength nor that of maximum growth rate are affected by the 
value of La*. However, the maximum growth rate itself is altered by the relative weight 
of inertial effects with respect to viscosity and capillarity. Therefore, the Laplace number 
is relevant when discussing time scales and growth rates, but not for critical or dominant 
wavelengths. 

The modes with A; > 1 correspond to the r > 0 region and are always stable as it 
is the case in the usual viscous lubrication approximation, but we want now to analyze 
whether there is any change in their behaviour when inertia effects are included. First, 
note that for each A: > 1, there is a single value of r > 0 (see Fig. 2a). This value 
of r could yield either ip = Q (blue line, upper branch) or \p\ < 7r/2 at the black line 
in Fig. 2b. Two different situations ensue. If r > r^, the solutions are on the p = Q 
(blue) line, i.e. the modes are monotonically damped, and two different values of q are 
admissible: one smaller and the other larger than go,b- At the point B = (rb,go,b)) both 
roots degenerate into a single one. For 0 < r < r;,, the roots are found along the black 
line, and the modes are oscillatory and damped. From Eq. (35), we find the wavenumber 
corresponding to point B as given by 

Thus, for 1 < A < there are two damped real modes, while for K > Kf, (r < r^) two 
oscillatory (complex) damped modes are possible with increasing frequency oscillations 
and stronger damping as K increases. 

In summary, the condition Q(r) = 0 (i.e., r real) yield three types of lines in the (go, p) 
plane, which can be classified as: 

1. ¥? = 0, which yields stable damped (real) modes, 

2. \p\ = 7r/2, which can be related to unstable purely growing (real) modes, and 

2>. p ^ 0,7r/2, that will produce stable oscillatory modes in time, i.e. complex conju¬ 
gate roots of H. 

The procedure to obtain the dispersion relation of the problem, i.e. Ll{K), for a fixed 
La* is as follows. Given a value K, we obtain the corresponding r (see Eq. (35) and 
Eig. 2a). Then, with this value of r, we find qq (e.g. using Eig. 2b). In the case of 
complex roots (black line) the corresponding value of is a consequence of requiring 
that 9(r) = 0 in Eq. (40). Once this is done, we obtain the full spectrum of modes as 
shown in Fig. 4. The dashed lines correspond to for the complex modes along the 
black line named C. 

We observe in Fig. 4 that La* strongly modifies some features of the complete dis¬ 
persion relation. For instance, it modifies the maximum, Ltmax, in the unstable region 
{K < I, > 0). Note that the product La*Llrnax grows with La* because ^max de¬ 
creases with La* with an exponent less than one (see Fig. 3b). Analogously, La* also 
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Figure 4: Real (solid lines) and imaginary (dashed lines) parts of 14 = multiplied 

by La* as a function of the wavenumber K for (a) La* = 1, and (b) La* = 10. 
The curves for 14^ > 0 and A' < 1 (unstable region) correspond to \ip\ = 7r/2, 
and those for 14^ < 0 and K < Ki, (stable region for damped modes) correspond 
to (/? = 0 


affects the minimum in the stable region with AT < 1. For A' > 1, La* only modifies the 
value of Kb (see Eq. (48)). 

In Fig. 5 we show a more detailed comparison of the dispersion curves for several 
La*’s, both on the real growth rates for unstable (14^ > 0) and stable (wr < 0) modes. 
Part a) shows that as La* increases the unstable modes have lower growth rates, but 
the wavenumber of the maximum growth is not altered, and remains at Kmax = 

For very small La*, the viscous dispersion relation is rapidly approached (see Eq (34)). 
Eigure 5b shows the stable region of the instability diagram (AT > 1). Eor 1 < A' < Kb), 
there two branches of modes that decay exponentially, a characteristic of the instability 
which is lost in the viscous approximation. Eor AT > 1, the viscous solution, Ll^is, is a 
fairly good approximation if AT < Kb, but fails for K around Kb- Clearly, this solution 
cannot describe the oscillating modes for K > Kb- 


4 Bidimensional flow: Linear stability analysis 

We consider here the full Navier-Stokes equation, Eq (1), in its dimensional form without 
the long wave approximation assumptions, i.e. the ratio e is not necessarily small now. 
Therefore, the small perturbations of the free surface are done on the velocity and 
pressure fields, and are expressed in terms of normal modes with a wavenumber k = 
(fca;,0,0) parallel to the substrate. Thus, we have 

5v = v(z) 

p = pQ + 5p = pQ+ pi{z) (49) 

h = ho + 5h = ho + (:e^^'^^^\ 
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Figure 5: Dispersion relations, Qrik), for some values of La*: (a) stable region, and (b) 
unstable region. The dashed line for La* = 0 is given by Eq. (34). 


where v = {u{z), 0, w{z)) and 5h is the Lagrangian displacement of the free surface. Note 
that, for small perturbations, we have ^ = w{l)/uj. Then, the Navier-Stokes equation 


at first order in the perturbations becomes 

p dt5v =—V 5p + fiA6v. (50) 

Since we assume incompressible flows, ik • v = —Dw, where D = d/dz. In order 
to reduce the number of variables, we eliminate the pressure terms, by taking the 2 ; 
component of the V x VxEq. (50). After some calculations, we obtain 

{D"^ - kl){D‘^ - k‘^)w = ^ (51) 

where k = kx, and 

= 1 + u]/{vk‘^), (52) 

or equivalently 

a; = (s^ — 1 ) vk"^. (53) 

The general solution of Eq (51) is 

w = Ai cosh(/i: 2 ;) + A 2 cosh(s kz) + A 3 sinh(/c 2 ;) + A 4 sinh(s kz), (54) 


where the constants Aj (i = 1, ...,4) are calculated by applying the following boundary 
conditions. 

First, we impose the no flow condition through the rigid substrate, 

w\z=o = 0- (55) 

Second, we shall assume that there is no slip at the substrate, = 0. Since the 

flow is incompressible, then 

Dw\z=o = 0 (56) 
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Third, the tangential stresses at the free surface should be zero, k • 5 • = 0, 

where S = —pi + pS/Sv + /x(V5v)^ is the stress tensor. By replacing here the perturbed 
quantities, Eq. (50), we hnd 

iD^ + k^)w\^^^^=0, (57) 

Finally, the normal stress at the free surface must satisfy the generalized Laplace 
pressure jump, 

=7C + n, (58) 

where C = —k'^Q is the first order curvature of the perturbed free surface. Since 


c 


w 

OJ 


z=ho 


(59) 


we have 

C - (60) 

h=ho 

Notice that the term in dU/dh plays a role that is analogous to that of pg in the Rayleigh- 
Taylor instability of a thin film. In order to obtain pi for this equation, we perform the 
scalar product of Eq. (50) by k, and using Eq. (53) we find 


, ,, dli 

{-pi + 2pDw)\^^. = — 


dh 


Pi 



(61) 


Then, by replacing this expression of pi ai z = ho into Eq. (60), we hnally write this 
boundary condition as 


1 

£2 


{k'^f 


D- 

CO 


, L»2 

= pi—-2- 

z=ho \ k^ 


Dw 


z=ho 


(62) 


where i is dehned by Eq. (20). 

From the above boundary conditions, Eqs. (55), (56), (57) and (62), it is possible to 
build up a matricial system to solve the four unknowns, Aj (i = 1,..., 4). Its determinant 
must be zero to avoid a trivial solution. This condition leads to 


[—4 (s + s^) + s (5 + 2s^ + s'^) cosh(iLe) cosh(s Ke) — 

(l + 6s^ + sinh(Are) sinh(s Ke)\ + 

La (AT^ — l) [s cosh(s Ke) sinh(iLe) — cosh(iLe) sinh(s Ke)] = 0, (63) 

with K and La dehned in Eqs. (25) and (31), respectively. 

This expression is the dispersion relation of the problem, since it implicitly gives s as a 
function of K. The values of co can be obtained through Eq. (53), which in dimensionless 
variables is 

/v 2 

Q = COT = „ (s^ — 1). (64) 

e-^La 

It can be shown that Eq. (63) is identical to that obtained in Kargupta et al. (2004) if 
slipping at the substrate is neglected once we take into account that a and f3 in their 
Eqs. (6) and (7) are Ke and sKe respectively (our s corresponds to their q). 
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In order to obtain the limit of Eq. (63) for e <C 1, note first that q in Eq. (29) of the 
long wave model is related to s in Eq. (52) by 

= (1 _ s‘^)K^e^, (65) 

and that Ke = 27r/io/A << 1 in this limit. In order to keep a meaningful value of q, 
|s| ;§> 1 is required, which means that q ~ iKse. Thus, with these ingredients in mind 
when analysing the limiting behavior of the dispersion relation given by Eq. (63), to the 
lowest meaningful order in e, we find 

La* K‘^{K^-l){ia.nq-q)=q^ (66) 

which is the same expression as given by the long wave model when Eqs. (35) and (36) 
are combined. 


5 Comparison between long wave and bidimensional 
models 


In this section we study the effects of La and e on the dispersion relation for the unstable 
region as given by the one dimensional (ID) long wave approximation and the bidimen¬ 
sional (2D) model. For the ID case, we focus on the solution of Eqs. (35) and (36) for 
(p = 7r/2, while for the 2D case we numerically solve Eq. (63) together with Eq. (64). 

In Fig. 6 we show the comparison between ID and 2D dispersion relations for given 
values of La (columns) and £ (rows). The inertial effects are shown along a given row 
(fixed e), with the first column being a viscous dominated flow, and the fourth column 
corresponding to inertia dominated cases. For small e, as in first row where e = 0.1, 
both dispersion relations are practically coincident for any value of La, as expected and 
shown analytically in Eq. (66). In general, the long wave model qualitatively predicts 
the same trends as the 2D model. However, for e as large as e = 0.5 (second row), 
the quantitative comparison certainly depends on La: the smaller La, the larger is the 
departure between both models, i.e. 2D effects become more important for flows with 
weak inertia. This effect is still more pronounced for larger e as seen in the third row 
for e = 1. Also note that, for fixed La, the position of the maximum shifts towards the 
left as £ increases. 

We focus now on the behavior of the maximum of the dispersion relations, since its 
analysis provides interesting insight on the effects of both inertia and aspect ratio. While 
the behavior for the ID model has been already described, the 2D model results can be 
obtained noticing that for D = Llmax'- 

^ dQ(s,K) (9D dQ ds 

dK dK dsdK ^ ’ 

Since the dispersion relation satisfies F{s,K) = 0, one can calculate 


ds 

dK 


dF 

dK 

dF 

ds 


( 68 ) 
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Figure 6: Dispersion relations, D as a function of K, for different values of La and e 
for the linearized problem: 2D model (N-S solution, blue dots) and ID long 
wave approximation with inertia (purple lines). For large La (strong inertia) 
we obtain similar results for both models. For small La (weak inertial effects) 
there are meaningful differences between both models (different e’s) 


Thus, using Eq. (64), it is possible to write 

which we shall not write in full for brevity. By solving this expression in conjunction 
with Eq. (63) we are able to obtain Llm and K^a as a function of both La and e. 

In Fig. 7 we show the wavenumber at the maximum growth rate, Km, as a function of 
La for several aspect ratios e’s, and vice-versa. Recall that for ID model, we simply have 
Km = l/\/2, independently of both La and e. For small La, the departure between both 
models can be very large if e is not very small. In fact, the value of Km can be reduced 
even up to 50% for e as large as e = 5 (see Fig. 7a). The difference remains also for 
large La, but it reduces for smaller e. This effect is clearly shown in Fig. 7b since, even 
if the departure increases for e increasing, it is smaller for larger La's. Therefore, the 
lubrication and the long wave approximations predict pretty larger distances between 
drops after breakup if the corresponding aspect ratio does not fulfill the requirement 
e <C 1. However, this discrepancy is smaller for larger La's. 

Figure 8 shows the maximum growth rate, Llm, as a function of La for several aspect 
ratios e’s, and vice-versa. The curves for ID model are obtained for the corresponding 
value of La* as given by Eq. (30). The difference in Llm between both models for small 
La can be very large if e is sufficiently large (see Fig. 8a). Instead, for large La both 
models agree in a power law decrease of Llm with differences that increase for larger e, as 
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Figure 7: Wavenumber at the maximum growth rate, as a function of: (a) La for 
different e’s, and (b) e for different La’s. The solid lines corresponds to 2D 
model, and the dashed line to the ID (long wave) model, = l/\/2 = 0.707. 


expected. The departures for small La and large e are also seen in Fig. 8b which shows 
how the ID curves with smaller La separate more and more from the corresponding 2D 
ones for smaller La’s. Thus, for small values of La the discrepancy in Llm between ID 
and 2D models can be very large even if e is not strictly much less than one (see also 
Fig. 10). 

6 Numerical simulations 

In order to analyse the validity range of the predictions of both LSA’s described above, 
we perform numerical simulations of the instability by solving the complete set of Navier- 
Stokes equations. Here, we use the two-phase flow, moving mesh interface of COMSOL 
Multiphysics. It solves the full incompressible Navier-Stokes equations using the Finite 
Element technique in a domain which deforms with the moving fluid interface by us¬ 
ing the Arbitrary Lagrangian-Eulerian (ALE) formulation. The interface displacement 
is smoothly propagated throughout the domain mesh using the Winslow smoothing al¬ 
gorithm. The main advantage of this technique compared to others such as the Level 
Set of Phase Field techniques is that the fluid interface is and remains sharp. The main 
drawback, on the other hand, is that the mesh connectivity must remain the same, which 
precludes the modelling of situations for which the topology might change. The default 
mesh used throughout is unstructured and has 2940 triangular elements (PI linear el¬ 
ements for both velocity and pressure). Automatic remeshing is enabled to allow the 
solution to proceed even for large domain deformation when the mesh becomes severely 
distorted. The mesh nodes are constrained to the plane of the boundary they belong to 
for all but the free surface. 

We adapt the same physical boundary conditions used above to the complete (nonlin- 
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(a) (b) 


Figure 8: Maximum growth rates, 0^, as a function of: (a) La for different e’s, and (b) 
e for different La’s. The solid lines corresponds to 2D model, and the dashed 
line to the ID (long wave) model (same colour implies same e in (a), and same 
La in (b)). The upper dotted line in both figures is the purely viscous {La = 0) 
growth rate, Llm^vis = 1/12 = 0.0833. In (a), the ID and 2D models for e = 0.1 
are graphically superimposed to the this value of Llm- 


ear) 2D problem. Thus, we write the kinematic condition as: 



(70) 


n being the external unit normal vector. Both surface tension and disjoining pressure 
exert normal stresses at the liquid-air interface 


5 • n = {aC — Kf{h)) n, 


(71) 


where C = — V* • n is the curvature of the free surface, = Ig • V is the surface gradient 
operator, and = I — nn the surface identity tensor. At the ends of the domain 
(x = 0 and x = d) periodic boundary conditions are applied for both the velocity 
field and shape of the free surface. On the liquid-solid interface, the the no-slip and 
no-penetration conditions (v = 0) are applied. 

Since we must have the same length scale in both x and z directions in the solution 
of the full non-linear N-S equations, we define now a slightly different dimensionless 
set of units than in LSA for the long wave approximation (see Eq. (25)). Thus, the 
dimensionless variables in the numerical simulations are given by: 


Z = Iz, X = ix, t = Tt, 


1 ~ 7 ~ 1 ~ 

u = —w, w = —w, p = —p, 
p pL t 


(72) 
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Figure 9: Time evolution of thickness profile for La 


1 and e 


1. We use Aq = 0.05. 


which yields the dimensionless form of Navier-Stokes equations 

dp f d‘^u d‘^u'\ 

dx dz^ j ’ 

dp / d'^w d‘^w\ 

(9z \ dx^ dz^ ) 

In particular, the boundary condition at the free surface, Eq. (71), becomes 


^ du du du\ 


La 


dw _ dw _ dw 


(73) 

(74) 


5 • n = (c- —f(h}) n. (75) 

\ go J 

In order to observe the evolution of the mode with maximum growth rate of the 2D 
model, we choose the length of the domain size in rr-direction as d = \m = . 

Note that this value is not coincident with (see Fig. 7). Thus, we use the following 
monochromatic initial perturbation of the free surface 

~ ( Sttx \ 

h{x,t = 0) = s + Aq sin ( j , (76) 

where Aq is a small amplitude (Aq = 0.05 in the present calculations). In Fig. 9 we show 
a time evolution of the thickness profile for La = 1 and e = 1 (we use (n, m) = (3, 2) 
and h* = 10“^ in all the following cases, unless otherwise stated). We carry on the 
simulation until the him becomes too close to /i*, where the numerical method is unable 
to converge, although continuation is sometimes possible by using automatic remeshing. 
We study the evolution of the instability by tracking the maximum and minimum 
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amplitudes of the free surface deformation by defining 


■^max (g 


max 1 
0<5<L 


h{i) 


min 1 
0<£<L 


h{i) 


(77) 

(78) 


These results are plotted in Fig. 10 for the same values of La used in Fig. 6, but e = 
0.5,1,2. The numerical non-linear solution of the problem shows that both Amax and 
Amin are practically coincident during a relatively long time of the evolution. Within 
the wide ranges of La and e shown in Fig. 10, this behavior is observed for at least two 
thirds of the total time required for the full development of the instability. This indicates 
that a linear models, like those presented previously, are relevant to describe the flow 
beyond the onset of the instability. 

In order to compare the numerical results with the linear models, we plot in Fig. 10 
the expected exponential behavior as, 

^ = 0.05e"(^“^°) (79) 

where a is given by the predicted growth rate for and To stands for a possible time 
shifting. For ID model, a is the corresponding growth rate for A^, i.e. a = 
which in general does not coincide with the maximum growth rate within this approx¬ 
imation. For 2D model, a = D^^(A^) = Llm, which is indeed the absolute maximum 
for this approach. Moreover, after separation of Amax and Amin, we expect that Amax 
remains closer to the exponential growth than Amin, which is more strongly affected by 
the presence of the substrate. This effect is certainly observed in the numerical results. 

Figure 10 shows that for small La, say La = 0.01 and 1, there is a very good agreement 
with the exponential behavior of the 2D model prediction (solid blue lines). In general, 
the ID model is not a good approximation, except for very small e, as expected. For both 
models, we use Tq = 0 since the behavior of Amax and Amin is of the exponential type 
from the very beginning. This type of growth is also observed for La = 100 and e = 0.5, 
but To / 0 is needed for large e, thus indicating the presence of a very early stage with 
slower (non exponential) growth. This effect is still more pronounced for La as large as 
La = 10^. In these cases, where there is still an acceptable agreement between the 2D 
model and the numerics for relatively large e. However, for this very large value of La, 
as £ is decreased neither ID nor 2D models are able to capture the actual evolution of 
the complete nonlinear problem. This issue deserves further investigation, which is out 
of the scope of the present paper and remains for future work. 

It is worth noting that the exponents n and m of the disjoining pressure in Eq. (3) do 
not a play a role in both linear analyses performed here. Their influence in this stage 
is somehow hidden in the length scale, £, defined in Eq. (20). However, some effects 
are expected in the numerical solution of the fully nonlinear N-S equations, since they 
appear in the boundary condition given by Eq. (75). Figure 11a shows a comparison of 
the time evolution of Amax and Amin for (n, m) = (3, 2) and (9,3), which are typical pairs 


Inertial and dimensional effects on the instability of a thin film • 2015 


page 19 of 25 









A. G. Gonzalez, J. A. Diez and M. Sellier 



T T T T 


Figure 10: Time lines of the amplitudes Amax (o) and Amin (a) with T = ie^ for different 
values of La and e. The lines correspond to the exponential behavior A = 
0.05 exp[0^(7" — To)], where Qm corresponds to the value given by either 
the 2D (solid line) or ID (dashed line) model. Tq = 0 except for the cases 
(Ta,e) = (10^, 0.5), (lO^l) and (10“^, 1). 


of exponents used in the literature (Schwartz, 1998). Clearly, both cases are practically 
coincident in the linear stage, and are in agreement with the linear 2D model. For larger 
times, the corresponding non-linear regimes strongly differ, thus leading to different 
breakup times, so that the effect of the exponents is limited to the short final non-linear 
stage. Similarly, Fig. 11b shows the same time lines for two different values of /i* and a 
given pair of (n, m). Also in this case, only the nonlinear stage of the evolution changes 
for different thicknesses /i*, without any significant change of the early linear stage. For 
fi* as small as = 10“^, no difference is observed neither in the linear nor in the 
nonlinear stages. This so because /i* becomes negligible with respect to ho. 

7 Summary and conclusions 

In this work, we have developed three different approaches to study the instability of 
a flat liquid thin film under partial wetting conditions, and subject to intermolecular 
forces (disjoining pressure): long wave ID model (with inertia), liner 2D model, and 
fully nonlinear numerical simulations. Firstly, we have extended the purely viscous 
analysis within the lubrication approximation to one where inertial effects are taken 
into account, which we call for brevity ID model. The LSA of this model shows that 
inertia does not lead to new regions of instability compared with the purely viscous case. 
Instead, it adds new stable modes: some which are exponentially decaying, and others 
which are damped oscillations. The former extend over the same range of the unstable 
modes and even beyond, while the latter appear for larger wavenumbers. In the unstable 
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(a) (b) 

Figure 11: Time lines of the amplitudes A^ax and Amin for La = 1 and e = 1 for: (a) two 
different pairs of exponents {n,m), (b) two different values of h^. Symbols 
indicate numerical simulations, lines the predictions of linear models. 


region of most interest here, we find that both the marginal wavenumber and that of the 
maximum growth rate do not change at all with the addition of inertia. However, the 
results clearly show that the growth rates of the instability decrease as inertial effects 
are stronger. The intensity of these effects is here quantified by a single parameter, 
namely the modified Laplace number. La* = Lae^. Therefore, the approximation can 
be applied only for large La, since e <C 1 is required for the approach to be valid. 

Secondly, we develop a LSA of the Navier-Stokes equation, so that the restriction 
of small aspect ratio, e, is no longer required. This calculation, called for brevity 2D 
model, is particularly useful to assess the accuracy of the ID model predictions. The 
main difference between these models is the way that inertia is treated. In linear ID 
model, the convective terms for the horizontal direction are still taken into account, 
while horizontal and vertical convective terms are neglected in the linear 2D model, but 
the viscous Laplacian term is now fully conserved for both directions. Thus, we have 
now two independent parameters to characterize the flow, namely La and e. The 2D 
model shows that the marginal wavenumber remains the same as ID model, and does 
not depend on La. However, unlike ID model, 2D model shows that Km is not constant, 
and decreases as La increases. This is an important result, since it shows that inertia 
can modify the distance between the final drops, which must be more separated with 
respect to the purely viscous case. 

With respect to the dependence of the growth rates with La, 2D model also shows 
that they decrease for increasing La, but the strength of the effect is greater than what 
is predicted by ID model. Interestingly, the discrepancies between both models decrease 
as La increases, i.e. for larger inertial effects. Note also that both models capture the 
main scaling of the dimensional growth rate, u), with the aspect ratio e. Thus, we can 
write, 

3 3 

0Ji{k) =—LLi{kt, La£^), uj 2 {k) =—kl 2 {k£; La,e) (80) 

T r 
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^ = h. / hg 


Figure 12: Function F{^) that determines the influence of h* on the characteristic length 
scale, £. 


where the subscripts 1 and 2 correspond to ID and 2D models, respectively. 

Finally, we are concerned now with obtaining a prediction of the both km and LOm as 
a function of the film thickness, /iq, for a given experimental configuration. In order to 
do so, we recall that (Israelachvili, 1992) 


K = 


A 


(81) 


where A is the Hamaker constant. Thus, the characteristic length, given by Eq. (20) 
can be written as 

e = F{0^, (82) 

where 

being ^ = /i*//io. The function F{^), which describes the effects of h*, is shown in Fig. 12 
for three usual values of (n, m). Interestingly, the cases with m = 3 and large n present 
a practically constant region for ^ < 0.5, which is a typical range in experiments. In 
these cases, we notice that £ oc h^, thus e oc hg ^ and La oc h^. Instead, if m < 3, say 
m = 2, the behaviour is different since F —)• 0 for decreasing /i*. For (3,2), we have 
i oc hg^^, thus e oc hg and La oc ■ 

These results shoud be taken into account when analyzing experimental data within 
a given hydrodynamic model. For Instance, the lubrication approximation would not 
become more valid as ho decreases (as it could be expected a priori) since e increases for 
thinner hlms. In fact, let us consider the data from the experiments with melted copper 
films on a Si02 substrate reported in (Gonzalez et ah, 2013). In this case, we have 
7 = 1.304:Kg/m^, /r = 0.00438Pas, and the experiments could be fitted with a purely 
viscous lubrication model using A = 2.58 10“^® J, h* = 0.1 nm and {n,m) = (3,2). 
Thus, we calculate the corresponding values of e and La for film thickness, ho, in the 
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Figure 13: Dependence of the dimensionless parameters £, La, La*, and the character¬ 
istic length scale, i, as a function of the film thickness, ho, for melted copper 
films 


interval (1,100) rem, as shown Fig. 13a. Note that even if inertial effects increase as ho 
increases, e decreases even faster, so that lubrication approximation assumptions apply 
for larger /iq’s (see also La* in Fig. 13b). Consistently, Fig. 13b indicates that the length 
i (proportional to the critical wavelength) increases with ho, so that wavelengths of some 
hundreds of nanometers should be expected for these film thicknesses. 



Figure 14: Wavelength of maximum growth rate. Am, and the corresponding growth, 
ujm, as a function of the film thickness, ho, for melted copper films 

In particular, we show the wavelength of maximum growth rate. Am, as well as the 
corresponding growth, ujm, as a function of ho in Fig. 14. The asymptotic power laws 
for large ho, given by the lubrication approximation where = l/\/2 and Dm = 1/12 
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are (see also Eq. (25)), 


A 


m — 





-S'- 

3/r hi 


-3 

0 • 


(84) 


These expressions are plotted as dotted lines in Fig. 14. Therefore, we conclude that 
both inertial and bidimensional effects are not significant if /iq ^ 20/i* , being pretty safe 
to use lubrication approximation results to describe the instability even for large La 
provided ho ^ h^:, as in the experiments reported by Gonzalez et al. (2013). However, 
for very thin nanometric films with ho < lO/i*, these effects should be taken into account, 
specially when analyzing the growth rates of the unstable modes. 
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